package com.numericalmethod.algoquant.model.ralph2009.calibration;

import java.io.BufferedReader;
import java.io.FileReader;
import java.util.Vector;

import com.numericalmethod.suanshu.matrix.doubles.ImmutableMatrix;
import com.numericalmethod.suanshu.matrix.doubles.matrixtype.dense.DenseMatrix;
import com.numericalmethod.suanshu.stats.descriptive.covariance.Covariance;
import com.numericalmethod.suanshu.stats.descriptive.moment.Variance;
import com.numericalmethod.suanshu.stats.timeseries.datastructure.multivariate.realtime.inttime.MultivariateIntTimeTimeSeries;
import com.numericalmethod.suanshu.stats.timeseries.datastructure.multivariate.realtime.inttime.MultivariateSimpleTimeSeries;
import com.numericalmethod.suanshu.stats.timeseries.linear.multivariate.stationaryprocess.arma.VARFitting;
import com.numericalmethod.suanshu.stats.timeseries.linear.multivariate.stationaryprocess.arma.VARMAXModel;
import com.numericalmethod.suanshu.vector.doubles.ImmutableVector;



public class Ralph2009Calibrator {
		
	String InputFileCSV = "data/stock_n_yield.csv"; // input CSV file for stock price and yield	
	
	Ralph2009CalibrationParameter[] cali_parm; 
	Vector<Double> vec1;
	Vector<Double> vec2;
    double[][] cov;
    double[][] timeLag1_coeff;
    double[]const_coeff;
    boolean bReadInputBefore=false;
	
	
	int iWindowSize;
	double x1[];
	double x2[];
	int iOffset;
	double del_t;
	
	double x0_cur;
		
	private void Init_All() {
					
		vec1 = new Vector<Double>();
		vec2 = new Vector<Double>();
						  
		x1 = new double[iWindowSize];
		x2 = new double[iWindowSize];
		
        cov = new double[2][2];
        timeLag1_coeff = new double[2][2];
        const_coeff = new double[2];		
                
        del_t = 1.0;  		
	}
	
	// with user input rolling window size 
	public Ralph2009Calibrator(int iInput_WindowSize) {
		
		iWindowSize = iInput_WindowSize;
	}
	
	Ralph2009Calibrator() {
		
		iWindowSize = 22; // without user input, we set it to 1 mo. default widow				
	}
	
	
	public void calibrate() {
		doEstimationWithSuanShu();
	}
	

	private void esimateNow(double [][]timeLag1_coeff, 
			double []const_coeff,
			double [][]cov,
			double delta_t,
			int iTimeOffset) {
		
		// first we compute the as cholesky decomposed
		// i.e. A = L*L' where L' is the conjugate tranpose
		//[s	0] X [s	 x1]   [s^2	  sx1]
		//[x1	x2]	 [ 0 x2] = [sx1	  x1^2+x2^2]
		// we can write the closed form for all the lower triangular parameters.
		
		cali_parm[iTimeOffset].sigma_s = Math.sqrt(cov[0][0]); // // Eq Appendix_param(6)  s^(0.5)
		cali_parm[iTimeOffset].sigma_x1 = cov[1][0]/cali_parm[iTimeOffset].sigma_s; // // Eq Appendix_param(7)   cov(1,0)/s
		cali_parm[iTimeOffset].sigma_x2 = Math.sqrt(cov[1][1] - cali_parm[iTimeOffset].sigma_x1*cali_parm[iTimeOffset].sigma_x1); // // Eq Appendix_param(8)   sqrt(cov(1,1)-x^2)
		
		// fitted VAR(1) becomes
		// Yt = {const} + {coeff}Y(t-1)
		// expand as
		// [Mt] = [(1-phi)u0*del_t]  + [-(1-phi)del_t+1 u1(1-phi)del_t] [M(t-1)]    Sigma*eps
		// [Xt]   [ 0             ]    [0               -alpha*del_t+1] [X(t-1)] + 
		//
		//   [a] + [c d]Y(t-1)
		//   [b]   [e f]
		double a,c,d,f;
		f = timeLag1_coeff[1][1]; // Eq Appendix_param(1)
		c = timeLag1_coeff[0][0]; // Eq Appendix_param(1)
		a = const_coeff[0];       // Eq Appendix_param(1)
		d = timeLag1_coeff[0][1]; // Eq Appendix_param(1)
		
		cali_parm[iTimeOffset].alpha = -(f-1.0)/delta_t;// Eq Appendix_param(3) 		
		cali_parm[iTimeOffset].phi = 1.0 - (1.0 - c)/delta_t; // // Eq Appendix_param(2)
		cali_parm[iTimeOffset].u0 = a/((1-cali_parm[iTimeOffset].phi)*delta_t);// Eq Appendix_param(4)
		cali_parm[iTimeOffset].u1 = d/((1-cali_parm[iTimeOffset].phi)*delta_t);// Eq Appendix_param(5)
		
		// no momentum calibration (i.e. mean reversion)
		cali_parm[iTimeOffset].u0_prime = a/delta_t; // Eq Appendix_param(11)
		cali_parm[iTimeOffset].u1_prime = d/delta_t; // Eq Appendix_param(12)
		
		cali_parm[iTimeOffset].alpha_prime = cali_parm[iTimeOffset].alpha; // // Eq Appendix_param(10)   same as mometum version						
	}
	
	
	private void doEstimationWithSuanShu() {				
		readInputAndSetup();		
		
		// do it from T-n time back to time 0.
		for (int ii=iOffset; ii >= 0; ii--) {
			createOneRollingWindow(ii);
			callSuanShuToSolve(x1,x2);
		    esimateNow(timeLag1_coeff, const_coeff, cov, del_t,ii);
		    print(ii);		    
		    
		    //CallSuanShuNelderMead(x1,x2,ii); // this is to call the Nelder Mead to estimate by finding the max likelihood
		}		
	}
	
	void readInputAndSetup() {		
		
		if (bReadInputBefore == false) {
		    readInputFile();
		}
		
		cali_parm = new Ralph2009CalibrationParameter[vec1.size()];
		
		for (int ii=0; ii < vec1.size(); ii++) {
			cali_parm[ii] = new Ralph2009CalibrationParameter();
		}

		
		if ((vec1.size() - 1) >= (2*iWindowSize-1)) {			
			iOffset = (vec1.size() - 1) - iWindowSize - (iWindowSize-1); // 1st starting pt for storing the fitting
		}
		else {
			System.out.println("Not enough data for window size= "+iWindowSize + "for data point = " + vec1.size());
			iWindowSize = vec1.size()/2; // (2*win_size-1) <= total_element-1
			System.out.println("Resized window  =" +iWindowSize);
			iOffset = (vec1.size() - 1) - iWindowSize - (iWindowSize-1); // 1st starting pt for storing the fitting
		}				
	}
	
	
	private double calM(Vector<Double> TmpVec1,
			  int iTmpOffset) {
		double rate_ret, M=0, discount;
		
		for (int ii=0; ii < iWindowSize; ii++) {
			rate_ret = (TmpVec1.elementAt(ii+iTmpOffset) - TmpVec1.elementAt(ii+iTmpOffset+1))/TmpVec1.elementAt(ii+iTmpOffset+1);
			if (ii > 0) {
				discount = Math.exp(-ii);
			}
			else {
				discount=1.0;
			}
			
			M += discount*rate_ret;
		}
		
		return M;
	}
	
	private void createOneRollingWindow(int iTmpOffset) {
		double u_avg=0.0;
		
		// we assume no of elem in vector is at lest 1 elem greater than windowSize
		// i.e. vec.size=5, windowsSize=4; because for a window of generating 4 return, we need 5 points.
		
		for (int ii=0; ii < iWindowSize; ii++) {
			u_avg += Math.log(vec2.elementAt(iTmpOffset+ii)); // sum up all first
		}
		u_avg = u_avg / iWindowSize;
		
		for (int ii=0; ii < iWindowSize; ii++) {
			//x1[ii] = vec1.elementAt(iOffset+ii);
			x1[ii] = calM(vec1,iTmpOffset+ii);
			x2[ii] = Math.log(vec2.elementAt(iTmpOffset+ii)) - u_avg ; // this is X (de-meaned yield), yeild - mean(yield)
			
			//System.out.format("x1[%d]=%f x2[%d]=%f %n",ii,x1[ii],ii,x2[ii]);
			if ((iTmpOffset ==0) && (ii== 0)) {
				x0_cur = x2[ii];
			}
		}				
	}



	private void callSuanShuToSolve(double []x1, 
			double []x2) {
		DenseMatrix matA = new DenseMatrix(iWindowSize,2);
		
		// Now, input from backward because SuanShu expect index 0 is time 0, and index 1 is time 1.
		int jj=0, ii;
		for (ii=iWindowSize-1; ii >= 0; ii--) {
			matA.set(ii+1, 1, x1[jj]);
			matA.set(ii+1, 2, x2[jj]);
			jj++;;
		}

		MultivariateIntTimeTimeSeries mts = new MultivariateSimpleTimeSeries(matA); 
		
        VARFitting var1 = new VARFitting(mts, 1); // fit as VAR(1)
                 
        VARMAXModel varmax1 = var1.getVARMAX();
        
        ImmutableVector d_mu = varmax1.mu(); // get the constant term
        ImmutableMatrix mat_phi = varmax1.AR(1); // get AR coeff        		
        
        double mu_arr[] = d_mu.toArray();
        
        const_coeff[0] = mu_arr[0];
        
        timeLag1_coeff[0][0] = mat_phi.get(1,1);
        timeLag1_coeff[0][1] = mat_phi.get(1,2);
        timeLag1_coeff[1][1] = mat_phi.get(2,2);
        
                
        Variance sig1 = new Variance(x1);
        Variance sig2 = new Variance(x2);
        Covariance sig12 = new Covariance(x1,x2); 
        
        cov[0][0] = sig1.standardDeviation()*sig1.standardDeviation();
        cov[1][1] = sig2.standardDeviation()*sig2.standardDeviation();
        
        cov[0][1] = sig12.correlation()*Math.sqrt(cov[0][0])*Math.sqrt(cov[1][1]);
        cov[1][0] = cov[0][1];
	}
	

	private void readInputFile() {
		
		Init_All();
		
		try {
			BufferedReader CSVFile = 
			        new BufferedReader(new FileReader(InputFileCSV));
		
			String dataRow = CSVFile.readLine(); // Skip first line.
			dataRow = CSVFile.readLine(); // Read first line.
			  
			double dTmpX;
			
		
			while (dataRow != null){
			    String[] dataArray = dataRow.split(",");
			   
			   for (int ii=1; ii < dataArray.length; ii++) { 
				   
				   dTmpX = Double.parseDouble(dataArray[ii]);
				   
				   if (ii == 1) { // col 1
					   vec1.addElement(dTmpX);
				   }
				   else 
					   if (ii == 2){ // col 2 
						   vec2.addElement(dTmpX);
					   }
			  //    System.out.print(dTmpX + "\t"); 
			   }
			//   System.out.println(); // Print the data line.
			   dataRow = CSVFile.readLine(); // Read next line of data.
			}		 
			CSVFile.close();	
		} catch (Exception e) {
		    System.out.println(e.toString());
		}		
	}
	
	

	private void print(int iTimeOffset) {
		System.out.println("time = "+iTimeOffset);
		System.out.println("SuanShu sigma_s = " + cali_parm[iTimeOffset].sigma_s);
		System.out.println("SuanShu sigma_x1 = " + cali_parm[iTimeOffset].sigma_x1);
		System.out.println("SuanShu sigma_x2 = " + cali_parm[iTimeOffset].sigma_x2);
		System.out.println("SuanShu VAR u0 = " + cali_parm[iTimeOffset].u0);
		System.out.println("SuanShu VAR u1 = " + cali_parm[iTimeOffset].u1);
		System.out.println("SuanShu VAR alpha = " + cali_parm[iTimeOffset].alpha);
		System.out.println("SuanShu VAR phi = " + cali_parm[iTimeOffset].phi);
		System.out.println("ut= "+(cali_parm[iTimeOffset].u0+x2[0]*cali_parm[iTimeOffset].u1));
		System.out.println();
	}
	
	
}
